I. Introduction 

Whether you want to plan an efficient bus route, optimize the positions of in- 
tegrated circuits on a silicon chip, or predict the lowest energy configuration of a 
cluster of atoms, the problem reduces to one of finding minima or maxima of spec- 
ified functions. The significance of the optimization field is reflected in the number 

ON ' of papers written on the subject. Over 300 articles on one method alone, simulated 

ON " 

annealing, have been published since 1988. 

There appears to be no one "perfect" algorithm that will solve every optimization 
problem. Instead, a host of complementary methods have evolved, each being suited 
to specific tasks [1]. For example, Brent's method or the golden section search can 



variable even if its derivative is unknown. Multidimensional functions are much more 
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be used if the function to be minimized (or maximized) has only one independent 
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difficult to optimize. If the first derivatives can be calculated, conjugate gradient 
methods or variable metric methods may be used. Otherwise, the downhill simplex 

D . 

method, direction-set methods, simulated annealing, or genetic annealing may be 

O 

used. None of these techniques can be universally implemented, and the method of 

>< 

5_i ■ choice varies with the details of the specific problem. 
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Simulated annealing is a particularly promising minimization technique [1,2]. It 
has, for example, proved effective in finding the global minimum of multidimensional 
functions having large numbers of local minima. As with other Monte Carlo based 
approaches, this method is well suited for implementation on both serial and parallel 
architectures. The mathematical problem, the optimization of a function, is "solved" 
by relating the task to a corresponding physical process, thermal annealing. As clas- 
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sical systems are slowly cooled, the configuration space density tends to condense in 
regions where the potential energy is small. If the cooling process is slow enough, 
then the system ultimately finds the physical arrangement that minimizes its po- 
tential energy Hence, by assigning the function to be minimized to be the analog 
of the potential energy and some control parameter as the analog of temperature, 
the global minimum of a function can be found by simulating the annealing process 
as the "temperature" is taken to zero. As in the case of the physical system, the 
rate of cooling is important in determining whether or not such annealing procedures 
ultimately find their way to the global minimum or become trapped in various local 
minima. Although this technique is relatively new, it has proved useful for a wide 
range of optimization problems [1]. 

Within the protein folding community, various quantum methods have been de- 
veloped in which the global minimum is found either by tunneling out of local minima 
or by removing uninteresting minima through smoothing [3]. While this work was 
in progress, Amara, Hsu, and Straub developed a minimization scheme for multidi- 
mensional potentials via an approximate solution of the imaginary time Schrodinger 
equation. The approximate wave function is comprised of a Hartree product of single 
wave packets. These packets are allowed to move, tunnel, expand, and contract in 
search of the global minimum. The quantum mechanics of the system is then relaxed. 
In the limit of h -> 0, the classical minimum if found [4]. 

We describe in Section II a new optimization approach, quantum annealing, that 
is closely related to its classical counterpart. Unlike the other quantum mechan- 
ical approaches discussed, this method does not require an approximation to the 
wavefunction. We illustrate its use with two examples in Section III. The first, a 
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pedagogical exercise, is designed to illustrate how quantum annealing avoids local 
minima. The second is a non-trivial determination of the lowest energy configura- 
tions of various homogeneous Lennard- Jones clusters. In Section IV we present our 
summary and discussion. 

II. Formal Development 

We begin by assuming that our task is the minimization of a specified, many- 
variable function. As in conventional annealing methods, it is convenient to view 
this function as being the potential energy of a hypothetical physical system, in 
contrast to the simulated annealing approach, however, it proves convenient to view 
our system as quantum-mechanical rather than classical in nature. 

Assuming that we can compute the average energy of our system as a function of 
its temperatue and quantum mechanical character, we can approach our objective, 
the minimum of the system's potential energy (point A in Fig. (1)), in a variety of 
ways. Starting at the arbitrary point D in Fig. (1), simulated annealing first turns off 
the system's quantum mechanics and then reduces the temperature of the resulting 
classical system to zero (path DBA in Fig. (1)). It proves useful, however, to invert 
the order of these two limits first reducing the system's temperature to zero and 
then following the resulting quantum ground state energy to its classical limit (path 
DCA in Fig. (1)). Since the fictitious mechanical system involved is an artificial 
construct, we are free to craft its character to suit our purposes. In particular, we 
can control the degree of its quantum mechanical behavior by varying the masses of 
its constituent particles. 



3 



The motivation for inverting the usual order of the limits is that, in a sense, it is 
easier to take the zero temperature limit of a quantum problem than a classical one. 
in particular, a variety of diffusion based, Monte Carlo methods are available that can 
be used to compute the ground state energy of general quantum-mechanical systems. 
Described in detail elsewhere [5], these methods are based on the observation that 
the Schrodinger equation written in imaginary time is isomorphic to the diffusion 
equation with a growth depletion term. In one-dimension the result is 
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where r = ith, E is a constant subtracted out for convenience and V(x) is the 
potential energy. 

The Diffusion monte carlo (DMC) method is one, relatively simple technique 
for treating such problems. It follows the evolution of a number of random walkers 
designed to move in such a way to simulate the diffusion, growth and decay processes 
in Eq. (1). Walkers in regions where E < V(x) are attenuated via first-order decay, 
while those in regions where E > V(x) undergo analogous growth. In typical 
applications, E is adjusted iteratively to maintain a steady state population. In 
terms of the Schrodinger eigenfunctions and energies, {0„}, and {E n }, respectively, 
the solution we seek is of the form 



y(x,T) = J2Mx)e- (En - Eo)T 



(2) 



The ground state wavefunction can thus be identified as the large r limit of ^(x,r) 
while the ground state energy is equal to E . One advantage of using DMC is that 
no knowledge of the wavefunction is required. Instead, it can be obtained from the 
final distribution of walkers in the DMC simulation. 

From Eq. (2), it is apparent that the rate of convergence of the DMC method is 
controlled by the gap between the ground and first excited state energies (AE). In 
fact, the time required for the DMC method to converge to the ground state is H/ AE. 
This is the same time scale needed for tunneling between two interacting potential 
minima. This suggests that DMC finds the ground state wavefunction through tun- 
neling. By gradually increasing the mass of the walkers thereby constraining the 
wave function, we can follow the ground state energy to its classical limit. 

III. Numerical Examples: 

To illustrate quantum annealing, we first consider the problem of finding the 
minimum of the one-dimensional function shown in Figures 2 and 3. In Figure 2, 
we observe the evolution of transient distributions of DMC diffusers at a fixed mass 
to demonstrate how the method utilizes tunneling to find the ground state wave 
function. Initially, all the diffusers are placed in the left well at a ground state 
energy equal to approximately half the height of the intervening barrier (line a). 
Early in the simulation, some diffusers tunnel through the barrier to the deeper, 
right well (line b). Later, the population of diffusers in the deeper well grows relative 
to the population in the left well (line c) until the energy converges to the ground 
state and an equilibrium distribution is reached (line d). When the ground state 
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energy is large compared to the difference bewteen the two wells, the distribution 
is almost symmetric. However since the right well is marginally deeper, the right 
population will be correspondingly larger (assuming the number of diffusers is high 
enough to resolve the slight difference). 

In Figure 3, we show the evolution of the equilibrium population of our model 
system as we turn off its quantum mechanical character (lines a-c) . For a small mass, 
quantum effects are large and the ground stae energy is relatively high compared 
with the slight energy difference between the two wells. In this case, the populations 
in both wells are essentially equal (line a). As we increase the mass, the quantum 
character is reduced, the associated ground state energy drops below the energy of the 
metastable minimum and the probability of being found anywhere other than near 
the global minimum quickly dwindles (line c). In actual practice, it is generally not 
necessary to completely converge to the ground state probability distribution before 
increasing the system mass. Once even a single DMC diffuser reaches a particular 
potential well, it has the ability to proliferate and, hence, sample that region. As 
with conventional simulated annealing, however, it is necessary to experiment with 
annealing rates to be sure that results converge to the global minimum independent 
of initial configuration chosen. 

In our second example, we consider the problem of determining the structure of 
N atom homogeneous Lennard- Jones clusters. The particles are assumed to interact 
pairwise via a two parameter interaction of the form eV(£), where £ is a dimensionless 
interparticle separation distance, r/a. The parameters e and a are the usual Lennard- 
Jones energy and length scale variables. Explicitly, the Lennard- Jones potential is, 



6 



v(o = 4(r 12 - r 6 ) 



(3) 



In terms of our reduced parameters, the quantum mechanical Hamiltonian for our 

system is, 

\ H = -\ ^ 2 Ev 2 + E v(Un) (4) 

n=l m<n=l 

where 77 is a dimensionless constant that controls the scale of the quantum mechanical 
effects, 

h 

For rare gases, i] ranges in magnitude from 0.425 for 4 He to 0.00995 for Xe. As 
rj — > 0, the ground state energy approaches its classical limit. 

Figure 4 shows that the calculated ground state energies for Lennard- Jones clus- 
ters are relatively simple functions of the parameter rj. It is not difficult to relate 
the limiting slopes and curvatures of the ground state energies in Fig. (4) to cor- 
responding zero point energies and anharmonic effects, respectively. The simplicity 
of the ground state energies is in marked contrast to the excited states where the 
77-dependence is relatively complex even for a 3 atom cluster [7]. We note in par- 
ticular that the ground state energy of the clusters varies smoothly with 77 through 
regions where the behavior of the clusters is changing from delocalized to localized 
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in character. Since rj appears in Eq. (4) as the coefficient of the highest order deriva- 
tive, it is a "singular perturbation," an indication that power series expansions in 
the parameter 77 may have a limited (possibly zero) radius of convergence. 

The results of the quantum annealing calculations for various cluster sizes are 
shown in Table I. The energies quoted were obtained by extrapolating a linear fit to 
the ground state energies for small (< 0.01) values of rj. The extrapolated potential 
minima found are in good agreement with the known minimum energy structures of 
the classical system [8]. We note that for all cases reported in Table I, the agreement 
becomes exact if we refine the geometries of the large mass DMC clusters using 
traditional methods, a relatively simple task once the quantum annealing methods 
have successfully located the vicinity of the global minimum. That we obtain the 
correct structures is non-trivial since the 13 atom cluster has 988 known minima, 
and the 19 atom cluster has on the order of 10 5 stable isomers [9]. We note that the 
minimum energies found using the above extrapolation method are upper bounds 
to the actual minima since we are using linear fits to approximate a curve that is 
concave downward (see Figure 4). Decreasing the interval used in the fits (rj max ) leads 
to a better extrapolation of the curve and hence a better estimate of the potential 
minima. For example, in the seven particle cluster, as r\ max is decreased from 0.05 
to 0.02 to 0.01, the estimate for the potential minimum for the goes from -16.469 to 
-16.495 to -16.505, respectively, the final value being that reported by Hoare and Pal 
[8]. 

IV. Discussion: 
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Any method of locating the global minimum must address the issue of local min- 
ima. Simulated annealing confronts this problem through the device of classical 
"thermal fluctuations." Quantum annealing and other methods [3,4] use derealiza- 
tion and tunneling to avoid metastable regions. By utilizing a quantum rather than a 
classical system, the present approach exploits a number of specialized ground state 
methods that are not available within classical problems. Quantum annealing has 
the further advantage of making knowledge of the wavefunction unnecessary. The 
physically different ways in which quantum and simulated annealing avoid local min- 
ima suggests that these type of approaches may complement each other in general 
optimization applications. 
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Abstract 

Quantum annealing is a new method for finding extrema of multidimensional 
functions. Based on an extension of classical, simulated annealing, this ap- 
proach appears robust with respect to avoiding local minima. Further, unlike 
some of its predecessors, it does not require an approximation to a wavefunc- 
tion. In this paper, we apply the technique to the problem of finding the 
lowest energy configurations of Lennard- Jones clusters of up to 19 particles 
(roughly 10 5 local minima). This early success suggests that this method 
may complement the widely implemented technique of simulated annealing. 
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